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Dilute granular flows are routinely described by coUisional kinetic theory, but dense flows require 
a fundamentally different approach, due to long-lasting, many-body contacts. In the case of silo 
drainage, many continuum models have been developed for the mean flow, but no realistic statistical 
theory is available. Here, we propose that particles undergo cooperative displacements in response 
to diffusing "spots" of free volume. The typical spot size is several particle diameters, so cages of 
nearest neighbors tend to remain intact over large distances. The spot hypothesis relates diffusion 
and cage-breaking to volume fluctuations and spatial velocity correlations, in agreement with new 
experimental data. It also predicts density waves caused by weak spot interactions. Spots enable 
fast, multiscale simulations of dense flows, in which a small, internal relaxation enforces packing 
constraints during spot-induced motion. In the continuum limit of the model, tracer diffusion is 
described by a new stochastic differential equation, where the drift velocity and diffusion tensor are 
coupled non-locally to the spot density. The same mathematical formalism may also find applications 
to glassy relaxation, as a compelling alternative to void (or hole) random walks. 



I. INTRODUCTION 

In spite of a venerable engineering literature Q, |^ |^ , 
the study of granular materials is attracting a growing 
community of physicists 'Bl ^ 0, |3| , who sense that ba- 
sic principles remain to be discovered |^, including "the 
chance to reinvent statistical mechanics in a new con- 
text" Of course, individual grains obey the laws of 
classical physics, but in many ways their collective behav- 
ior defies standard assumptions in statistical thermody- 
namics and hydrodynamics. Static properties of random 
packings, such as geometrical correlation functions [TH 
0, the jamming transition fTsL ll^ fl^ ll^ fl^ ll^ . and 
local force distributions 0, l20ll2lLl22ll23ll24ll25| are far 
from fully understood, but granular flows may pose even 
more fundamental open questions. 

Naturally, the most attention has focused on the 
regime of fast, dilute flow because it is closest to the 
familiar molecular-fluid state. In this regime, grains to 
undergo simple random walks due to binary (or possibly 
many-body) collisions, which differ from those in normal 
fluids only by being inelastic. From these postulates, 
continuum equations for hydrodynamics and heat trans- 
fer can be formally derived from various modiflcations 
of Boltzmann's kinetic theory of gases I27L I2M I29I 
|3£, ^1, 32, 33] . Due to inelastic collisions, external forc- 
ing must supply enough kinetic energy to approximately 
satisfy the assumption of thermal equilibrium. The re- 
sulting deflnition of a time-dependent, local "granular 
temperature" in terms of velocity fluctuations is some- 
what controversial, since classical statistical mechanics 
may break down with inelastic collisions [lO| . For exam- 
ple, non-extensive Tsallis statistics js^ can be derived 
from the simple "inelastic Maxwell model" of a granular 
gas [3^. In any case, it is clear that kinetic theory breaks 



'Electronic address: "baza nt@math.mit.edul 



down with increasing density [3^ l3^ , when long-lasting 
many-body contacts are formed. Due to gravity, such is 
usually the case in granular materials. 

In seeking new and relevant physics, therefore, it seems 
more fruitful to focus on the opposite regime of slow, 
dense flow — closer to the "granular solid" than the 
"granular gas" y. Mean flow flelds in draining silos, 
incline chutes, and shear cells have been studied exten- 
sively, but a simple unifled description remains elusive. 
Classical continuum models of dense flows are based on 
plasticity theory from soil mechanics 0, Although 
they remain popular in engineering, these models predict 
complicated patterns of velocity discontinuities ( "rupture 
surfaces" ) , consistent with some experiments [s^ |3^ 0| 
but not others IstI lill li^. liM lil ] , and can lead to vi- 
olent instabilities [4^ Hg . On the other hand, another 
class of early models for silo drainage (discussed below 
in section) involves o nly simple geometrical considera- 
tions m lHlil m Ha Recent attempts to de- 
velop continuum models for dense flows have introduced 
a remarkable variety of physical postulates, such as a 
temperature-dependent viscosity |52j |. density-dependent 
viscosity [531 Im ) , non-local stress propagation along hy- 
pothetical arches [s^, self-activated shear events due to 
non-local stress fluctuations [s^ |53, coexisting liquid 
and solid micro-phases governed by a Landau-like order 
parameter [sM l59j| . and shear- induced local rearrange- 
ments mediated by free- volume kinetics 60, 6IJ. 

Putting aside the question of which of these contin- 
uum models best describes the mean velocity, we focus 
here on the stochastic motion of individual particles out- 
side the coUisional regime, for which no theory is avail- 
able. Even in experiments, little attention has been paid 
to this fundamental issue, perhaps due to the practical 
difflculty of tracking individual particles in dense flows. 
Due to recent advances in digital video imaging, however, 
thousands of particles near the wall of a dense granu- 
lar flow can now be tracked simultaneously with very 
high resolution [s^, US (well beyond previous stan- 
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dards llll |^). Combining these techniques with 
magnetic resonance imaging and x-ray tomography can 
also yield three-dimensional information |65]. Complete 
microscopic motion in the bulk can also be obtained 
from granular dynamics simulations, now involving up 
to hundreds of thousands of particles in three dimen- 
sions m m mil m These revolutionary ca- 
pabilities should aid in developing a statistical theory of 
dense granular flows. 

In some sense, work along these lines has already be- 
gun with attempts to define a meaningful temperature for 
dense granular systems. The notion of "Edwards temper- 
ature" rZlIlTa . [73 , analogous to structural temperatures 
in glasses Q , has attracted increasing attention over the 
past decade. The Edwards temperature differs from its 
Boltzmann counterpart in that the entropy is based on 
static random packings at a given volume, rather than 
thermal configurations at a given energy. As such, it 
may be related to the compaction of powders under vi- 
bration 0,123 ■ Some support for the idea has come from 
recent molecular dynamics simulations verifying an effec- 
tive Einstein relation between mobility and diffusivity in 
an artificial shear flow |63|- Nevertheless, it is not clear 
how the Edwards equilibrium ensemble relates to hydro- 
dynamics and diffusion (if at all) since no analog of the 
Boltzmann equation of kinetic theory has been proposed. 
This would seem to require replacing the usual hypoth- 
esis of a random walk, due to ballistic motion between 
instantaneous collisions, with a new dynamical mecha- 
nism. 

Thus we arrive at the central question of this work: 

Is there a simple analog of the random walk 
for a dense granular flow? 

Roughly a century after random-walk theory was first 
developed, the original idea of a single entity, such as 
a mosquito , sound-wave amplitude [t^, IZg] , Brown- 
ian particle '79', 'so] , or stock price '81] , undergoing ran- 
dom displacements remains in most subsequent general- 
izations [s^. IsM l84l |. For example, the "continuous-time 
random walk" of MontroU and Weiss ^851, '86'| and re- 
lated models 83.1, which describe anomalous diffusion in 
semi-conductors [s^j , turbulent flows js^l and many other 
systems 89] via a random waiting time between steps, 
also assume an independent random walker, unaffected 
by its neighbors. The same is true of the "persistent 
random walk" of Fiirth |93 and Taylor [93| which intro- 
duces auto-correlations between steps, e.g. to model a 
transition from short-time ballistic to long-time diffusive 
motion 92, 93]. 

Of course, there are correlations between different ran- 
dom walkers in any condensed-matter system (or finan- 
cial market 113]), but in most cases it is reasonable to 
assume that the diffusion of a single particle is indepen- 
dent of all the others. This universal independence of 
single-particle trajectories is often related to attaining 
the thermodynamic limit, in which an enormous num- 
ber {N w 10^'^) of particles undergo very frequent colli- 



sions (w = lO"*^^ Hz) at finite temperature (in the classical 
sense). As a result, short-range many-body interactions 
have little affect on tracer diffusion at macroscopic length 
and time scales, aside from influencing parameters, such 
as the diffusivity. 

In stark contrast, the "cages" of neighboring particles 
in a dense granular material may remain at least par- 
tially intact for the entire duration of slow flow 37], due 
to strong internal energy dissipation and much smaller 
system sizes (e.g. N ^ 10^). Even at macroscopic length 
and time scales, therefore, particles must diffuse together 
in a cooperative fashion, and it is not clear what kind of 
statistics might describe their collective motion. Boltz- 
mann's kinetic theory of gases does not appear to be the 
appropriate starting point. 

In this article, we propose a simple mathematical 
model for cooperative diffusion and apply it to slow gran- 
ular drainage from a silo, as shown in Fig. ^ After sum- 
marizing our results in section we begin by showing 
that the only existing model for granular diffusion, based 
on the concept of diffusing voids, is fundamentally flawed 
and draw clues from experiments as to what could be 
missing. In section Hvl we introduce the general concept 
of a "spot" of cooperative diffusion, and in section IVI we 
show that it is suffices to predict a wide variety of ex- 
perimental data. In section IVll we briefly discuss Monte 
Carlo simulations with the model and spot-based multi- 
scale algorithms. In section IVlII we derive general con- 
tinuum equations for tracer diffusion in the spot model, 
focusing in section lVlIII on the case of granular drainage. 
We conclude in sections IIXIIXII by discussing how the 
theory might have broad applicability to cooperative dif- 
fusion in amorphous materials. 



II. SUMMARY OF RESULTS 

As a guide for the reader, we list our main results: 

1. A new microscopic picture of dense granular flows 
is presented in Figs. IHI and |7| Rising spots of free 
volume cause strongly correlated random displace- 
ments of blocks of neighboring particles. Each ran- 
dom spot path corresponds to a thick chain of rep- 
tating particles. 

2. In section IIVI simple calculations show how the 
spot hypothesis resolves the "paradox of granular 
diffusion" illustrated in Fig.Q] Particles diffuse sev- 
eral orders of magnitude more slowly than free vol- 
ume. 

3. In section^ predictions of the Peclet number, dif- 
fusion length, and cage-breaking distance, consis- 
tent with experimental data, are made by assum- 
ing volume fluctuations of order one percent. The 
vertical diffusion length of free volume, which is dif- 
ficult to measure experimentally, is also predicted. 



FIG. 1: Sequence of experimental photographs of a d = 3mm glass beads draining from a quasi-two-dimensional silo (courtesy 
of A. Samadani and A. KudroUi jijl '). Particles are colored according to their initial positions in equally spaced horizontal 
layers when the flow began. 



4. Direct experimental evidence for spots, through 
spatial velocity correlations, is shown in Fig. |S1 

5. The necessity of weak repulsive interactions be- 
tween spots naturally explains experimental den- 
sity waves, sketched in Fig.|51 which propagate up- 
ward in wide funnels and possibly downward in nar- 
row funnels. 

6. Monte Carlo simulations with the simplest version 
of the model, shown in Figs. ^| and ^2 accu- 
rately describe tracer-particle dynamics and some 
aspects of many-body correlations, although un- 
physical density fluctuations grow in regions of high 
shear. 

7. The density may be stabilized by minor relaxations 
within spots to enforce packing constraints, which 
amounts to a very efficient multiscale simulation 
algorithm for granular flow. These extra fluctua- 
tions, sketched in Fig. ^1 may be responsible for 
the nontrivial (sub-ballistic) super-diffusion seen in 
experiments at scales less than a particle diameter. 

8. In section IVIII the mean-field continuum limit 
of the model is analyzed, starting from non-local 
stochastic differential equation, Eq. (^51) . The 
result is a Fokker-Planck equation for tracer diffu- 
sion, Eq. (|30|l . whose coefficients in Eas. (|33|) and 
H37|l are coupled non-locally to the spot (or free- 
volume) density. 



9. A general tensor relation between transport coeffi- 
cients, Eqs. (|4(JII . is derived, which relates spatial 
velocity correlations and the diffusivities of parti- 
cles and free volume. 



III. TOWARDS A STATISTICAL THEORY 
A. The Kinematic Model for the Mean Velocity 

In the 1970s, engineers developed a very simple "Kine- 
matic Model" for the mean velocity field in a steadily 
draining silo Ij, such as the flow illustrated in Fig.Q] It 
is based on the foUowingconstitutive law, suggested by 
Nedderman and Tuzun |4!^ . 

u = bV±v (1) 

which postulates that the horizontal velocity, u, is pro- 
portional to the horizontal gradient, Vj^, of the down- 
ward velocity component, v (i.e. the shear rate). We re- 
fer to the constant of proportionality, b, as the "diffusion 
length" of free volume for reasons soon to become clear. 
The intuition behind this relation is that particles should 
tend to drift horizontally toward regions of faster down- 
ward flow, which are less dense and more accommodating 
of newcomers. Assuming that these density fluctuations 
are small enough to justify steady-state incompressibility 
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yields an equation for the downward velocity, 

which is simply the diffusion equation, where the vertical 
coordinate, z, plays the role of time. As such, "initial 
conditions" must be specified for the downward velocity 
at the bottom of the silo, and information propagates 
upward. 

For example, a point source of velocity (a narrow ori- 
fice) at z = in a quasi-two dimensional silo produces 
the basic solution, 



v{x, z) = 



y/Anbz 



(3) 



for the case of constant b. The flow is effectively confined 
to a region of parabolic streamlines due to the diffusive 
scaling. Ax cx -y/z. As show in Fig. ^ this prediction is in 
rather good agreement with experimental measurements 
of the bulk flow (well below the free surface) for large, dry 
grains 

Ei|4ll|4a,IM|4d|- (For fine powders and soils, the 
model breaks down, presumably due to cohesion between 
grains and/or the influence of the interstitial fluid 0,^3-) 

The fact that a single fitting parameter, b, suffices to 
reproduce the entire flow field reasonably well should be 
viewed as a major success of the Kinematic Model. Nev- 
ertheless, the model fell from favor in the 1980s, in part 
because the diffusion length seems hard to predict a pri- 
ori. Although b is consistently on the order of a several 
particle diameters, its precise value depends somewhat 
on the particle size and flow symmetry |l| . (It also varies 
with the humidity in cases outside the pure "granular" 
regime ^2-) ^ Gaussian flt to the downward velocity 
profile also seems to require a larger value of b in the up- 
per region of the silo, where the flow is more plug-like. 
This was also seen as problematic, although the model 
does not really require that 6 be a constant. 

Perhaps a deeper reason for the dissatisfaction with the 
Kinematic Model is the natural prejudice against a the- 
ory seemingly devoid of "physics" . The model depends 
only upon geometry, via the diffusion length and the silo 
shape, and not on any of the usual physical variables, 
such as momentum, mass, energy, temperature, stress, 
etc. In contrast, more complicated models from soil me- 
chanics involve a stress-based yield criterion, which gen- 
erally leads to bulk discontinuities in stress and velocity 
("rupture zones") iXifiS]- It is not clear that these mod- 
els provide a better overall description of silo drainage, 
but in any case it seems worth taking a fresh look at 
the simpler Kinematic Model and its possible microscopic 
justification. 



B. The Void Model for the Mean Velocity 

Although Nedderman and Tuzun introduced the 
macroscopic view of Equation Q as a constitutive rela- 
tion and Eq. |2Jl as conservation law 43] , both equations 




FIG. 2: The Void Hypothesis I. A single particle drops in- 
dependently from one available cage to another J48|,J^JM 
by (equivalently) exchanging with a rising void |47l Isil l97l . 

As a result, cage breaking occurs at the scale of the particle 
diameter. 
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FIG. 3: The Void Hypothesis II. On a lattice jH, a void 
trajectory emanating from the orifice causes a reptating chain 
of downward particle displacements (thick lines) , each at the 
scale of the particle diameter. 



had been derived years earlier by Litswiniszyn | 48ll49ll5(]j 
and Mullins 1^ from a statistical model which re- 
mains, to this day, the only microscopic description of silo 
drainage. The basic hypothesis is sketched in Figure El 
Litwiniszyn first suggested in 1958 that individual hard- 
sphere particles diffuse downward through a fixed array 
of available "cages" , which yields the Kinematic Model 
for the mean velocity in the continuum limit 48, 49, 50l|. 

Independently, Mullins in 1972 derived the Kinetic 
Model from the hypothesis that that particles move pas- 
sively downward in response to the upward diffusion of 
"voids" emanating from the silo orifice 01 Al- 
though formally equivalent to Litswiniszyn's particle- 
based model [SJI, Mullins' void-based model introduces 
a fundamentally different perspective, analogous to va- 
cancy diffusion in crystalline solids, which we will find 
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quite suggestive. (Note that the same mechanism of "hole 
diffusion" was first proposed in the free- volume theory of 
molecular liquids and glasses [sE HBl , discussed below in 
section 1x1') 

The Kinematic Model is easily derived from the contin- 
uum limit of the Void Model, assuming that voids do not 
interact. According to the Central Limit Theorem [83l| . 
the void concentration, p^,, (or probability density) tends 
toward a Gaussian profile satisfying the diffusion equa- 
tion, 



oz 



(4) 



where the parameter, &, is well-defined microscopic quan- 
tity, the void diffusion length (usually assumed to be con- 
stant). The mean particle velocity. 



Pp (u, - 



-Wo (-6 V±/9„,pi,), 



(5) 



is obtained by equating the particle flux density (on the 
left) with minus the void flux density (on the right), 
where wo is the mean upward drift velocity of the voids 
(related to the flow rate), and Pp is the particle density, 
which is roughly constant in a dense flow. Equations 
and 101 follow from Equations 0) and Q. 

When physicists, Hong and Caram, revisited the Void 
Model in 199f , they performed computer simulations of 
voids undergoing directed random walks on a lattice , 
as illustrated in Fig. |31 Although they did not perform 
any mathematical analysis, the void diffusion length is 
easily calculated for any regular lattice: 



h = 



Var(Ax„) 
2dhAz„ 



(6) 



where Ax^ is the random horizontal displacement (in 
dh = 2 horizontal dimensions) when a void moves up by 
Azy. For example, assuming isotropic transition proba- 
bilities to nearest-neighbor sites in the next lattice plane 
upward, we find b = for the face-centered cubic 

and hexagon close-packed lattices and b = (i/4\/2 for the 
body-centered cubic lattice. More generally, for any con- 
ceivable lattice approximating the random close packing 
of hard spheres, the Void Model predicts 6 ^ d, in con- 
trast to experimental measurements which always yield 
b > d, e.g. b = 1.3d_p3: b = 2.3d HI, 6 = 3.5d 
H3 1 and b = 2d- id |4l|. This suggests that the micro- 
scopic picture of the Void Model is somehow flawed, even 
though the mean velocity profile is quite reasonable. 



C. Particle Dynamics in the Void Model 

More serious problems with the Void Model become 
apparent when one considers diffusion and mixing, which 
it seems has not previously been done. MuUins and 
Litwiniszyn were content to use the Void Model as simply 
a means to derive the continuum equation Likewise, 




FIG. 4; Simulation of the experiment in Fig. fusing the Void 
Model on a triangular lattice, which results in far too much 
mixing. (Courtesy of Chris Rycroft.) 



in spite of doing discrete simulations. Hong and Caram 
only measured the mean velocity profile. 

To quantify the diffusion of particles, we consider ge- 
ometrical single-particle propagator, Vp{'x.\z,xo, zq), the 
conditional probability density of finding a particle at 
horizontal position, x, after it has fallen to a vertical po- 
sition, z, from an initial position, (xojZq)- For a steady 
mean concentration of voids on a lattice, p„(x, z), it is 
straightforward to take the continuum limit of the exact 
lattice dynamics (as will be described elsewhere) to ob- 
tain a partial differential equation for 7'p(x|z, Xq, zq) at 
length scales much larger than the particle diameter: 



dVp 

dz 



^bV^-{Vi_Vp 



VpVj 



^OgPv 



(7) 



This equation may also be derived heuristically as fol- 
lows. In a uniform flow [p^ = constant), a particle does a 
directed random walk downward with precisely the same 
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diffusion length as the voids going upward. As in Eq. Q 
for the void concentration, this yields the first two terms 
in Eq. lO, i.e. the diffusion equation with "— z" acting 
like "time" . In a nonuniform flow, the third term rep- 
resents a horizontal bias of the particle random walk to 
climb the horizontal gradient of the void concentration, 
which ensures that the propagator is eventually localized 
on the source of voids (the silo orifice). The "advection 
velocity" , 6Vj_ logpv, is given by the ratio of (minus) the 
horizontal void flux to vertical void flux, ±py/ py, ac- 
cording to the right hand side of Eq. The task of 
characterizing solutions to the new equation 10 is left 
for future work. 



D. The Paradox of Granular Diffusion 

Here, it suffices to note that the particle dynamics in 
the Void Model is firmly contradicted by experiments. In 
a uniform flow, Equation (|7|) predicts that particles fol- 
low statistically identical geometrical trajectories while 
falling as voids do while rising, with the same diffusion 
length, b. Figure ^ provides an elegant visual demon- 
stration that this prediction is incorrect. The Void Model 
would predict that a particle starting near the top of 
the silo would diffuse horizontally within an inverted 
parabolic region of the same curvature as the parabolic 
flow region (set by the diffusion of voids). Therefore, a 
collection of particles starting near the top of the silo 
inter-diffuses and mixes across most of the flow region 
before falling to the middle of the silo, as shown in 
Fig. 21 Instead, experimental particles mix much less 
while draining, even in the non-uniform flow region, as 
the initially sharp interfaces between layers of different 
colors remain relatively intact during drainage. 

The failure of the Void Model has a clear microscopic 
origin. Particle-tracking experiments on silo drainage 
(motivated by the theory below) have recently shown 
that the "cage" of neighbors of a particle is preserved 
over surprisingly long distances, on the order of hun- 
dreds of par ticle diameters, even comparable to the sys- 
tem size |37| . This experimental result violates the cen- 
tral hypothesis of the Void Model, illustrated in Fig.|21 in 
any of its forms. On a lattice [o^l, when a particle moves 
by interchanging with a diffusing void, it loses roughly 
half of its neighbors. As this process is repeated, the 
particle is likely to lose all of its original neighbors after 
falling by only a few times its own diameter. This quickly 
results in a rather unphysical degree of mixing, as shown 
in Figs. Eland El 

This same conclusion also applies to the early engineer- 
ing theories, which are more vague about the microscopic 
dynamics 1]. Litwiniszyn AQ, 5QJ postulated that 
particles jump downward from one fixed cage to another, 
but this implies that a particle completely breaks free of 
its neighbors with each random displacement. Similarly, 
MuUins 1^ 1^ postulated that individual particles 
move by exchanges with particle-sized voids, indepen- 



dently from the random motions of neighboring particles, 
and this again implies a cage-breaking length compara- 
ble to the particle diameter. In reality, particles diffuse 
in a cooperative fashion, somehow managing to preserve 
their cages of first neighbors over rather large distances. 

How is it possible that the Void Model is able to de- 
scribe the macroscopic flow while so poorly describing 
the microscopic dynamics? This illustrates the danger 
of judging a statistical hypothesis based on only mean 
quantities, as some authors have argued that the success 
of the Kinematic Model lends support to the Void Model. 
The "paradox" of dense granular flow is that the diffu- 
sivity of particles is several orders of magnitude smaller 
than the diffusivity of free volume. 



E. More Clues from Experiments 

A closer look at experimental images suggests the key 
to resolving this paradox. As shown in fig.Oa), an inter- 
face between differently colored (but otherwise identical) 
particles deforms considerably in a non-uniform, dense 
granular flow, and yet individual particles rarely (if ever) 
penetrate from one region into the other. This is in stark 
contrast to diffusion in normal liquids and gases, where 
the independent random walks of tiny molecules cause 
an initially sharp interface to become smoothly blurred 
over time, even in the absence of a flow. Instead, we see 
that the granular interface in the center panel of Fig. [31 
develops some waviness at the scale of several particle 
diameters. This suggests that particles somehow diffuse 
cooperatively in cage-like blocks, mostly staying together 
with their nearest neighbors. 

This intuitive notion is consistent with other experi- 
mental indications of an important length scale of several 
particle diameters, which is missing in the Void Model. It 
is well known that a silo will not begin to drain until the 
orifice is at least several particles wide, and mechanical 
blocking (due to arching) can occur at somewhat larger 
orifice widths . The empirical Beverloo correlation [o^ 
implies that the extrapolated outlet diameter where the 
flow-rate vanishes is roughly 1.5 particle diameters, which 
has been explained in terms of a controversial empty an- 
nulus near the edge of the orifice 0. For our purposes, 
the salient point is that a granular material does not 
drain one particle at a time, so it is impossible to inject 
individual voids at the orifice. Since packing constraints 
must be enhanced in the bulk compared to the orifice, 
it seems highly improbable that voids could form and 
propagate in the interior of the silo. 

On the other hand, if particles tend to preserve their 
cages while diffusing in the bulk, then they should only 
pass through the orifice in correlated groups. Particle- 
tracking measurements of silo drainage show some signs 
fiow intermittency until the opening is at least six parti- 
cles wide [sJl, consistent with the point where the Bev- 
erloo correlation becomes acceptable T] . This again sug- 
gests that extended blocks of particles must be allowed 
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FIG. 5: Left: Experimental photograph of the initial condition of a quasi-two-dimensional silo with a vertical interface between 
black and white glass beads (by Jaehyuk Choi [s^ *). Center; The interface (in the left half of the silo) at a later time in the 
experiment. Right: A simulation of the same situation using the Void Model on a square lattice (by Guaqueta [q^*). Length 
is measured in particle diameters (d = 3mm). 



leave the silo at the same time. Thinking in terms of 
the Void Model, it could also mean that some extended 
entity causing motion is released into the bulk, once the 
orifice reaches a critical width. 

This idea may seem strange, but it has a chance of 
being correct. The Void Model does reproduce one cru- 
cial feature of particle dynamics in slow, dense drainage 
experiments p37| : All fluctuations are independent of the 
flow rate, aside from a simple rescaling of time. One 
implication is that the Peclet number. 



Pe 



Ud 

Id' 



(8) 



is independent of the flow rate, or equivalently the diffu- 
sivity, D, is proportional to the local mean velocity, U. 
This strongly suggests that diffusion and advection are 
caused by the same physical mechanism, as in the Void 
Model. The common mechanism, however, cannot be the 
propagation of a void or a cage-breaking event, which im- 
plies Pe « 1. Instead, it seems to involve the cooperative 
motion of neighboring particles, which somehow results 
in much less diffusion, Pe 3> 1. 



IV. THE SPOT HYPOTHESIS 

A. A Mechanism for Cooperative Diffusion 

In developing a microscopic theory, we are constrained 
by the fact that the mean velocity is fairly well described 



by the Kinematic Model, which seems to imply the diffu- 
sion of some entity causing motion. Let us suppose that 
such an entity exists, but, since it cannot be a void, we 
will call it a "spot". The kinematic parameter is then 
approximately set by the spot diffusion length. 



Var(Ax5) 
2dhAz„ 



(9) 



where Ax^ is the random horizontal displacement of a 
spot as it rises by Az^- 

What exactly is a spot? Since spots are injected at the 
silo orifice as particles drain out, it is natural to assume 
that each carries a certain amount of free volume. Rather 
than being concentrated in a void, however, this volume 
should correspond to a slight excess of interstitial space 
spread across an extended region, typically larger than 
a particle, as shown in Fig. We expect the size of 
a spot to be at least three particle diameters since it 
should induce the cooperative motion of a particle and 
its nearest neighbors. 

The simplest dynamical hypothesis, illustrated in 
Fig. 13 is that a spot causes all affected particles to move 
as a block with the same small displacement in the oppo- 
site direction. As the spot follows its random trajectory 
upward, a thick chain of particles cooperatively reptates 
downward by a very small distance, much less than the 
particle diameter. The subsequent passage of other spots 
causes each particle to randomly reptate in many differ- 
ent thick chains, as shown in Fig. 
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FIG. 6: The Spot Hypothesis I: Locahzed spots free inter- 
stitial volume, extending several particle diameters, diffuse 
upward from the silo orifice as particles drain out, causing 
the reverse reptation of thick chains of particles. 



of dense drainage. 

This microscopic picture is fundamentally different 
from previous theories, which assume that particles un- 
dergo independent random walks. According to Kinetic 
Theory, a particle moves ballistically between instanta- 
neous, randomizing collisions, and in the Void Model it 
jumps from cage to cage. In both cases, the particle 
loses a most of its neighbors in a single displacement, 
and thorough mixing occurs at the scale of several parti- 
cle diameters. 



B. Correlations Reduce Diffusion 

A back-of-thc-envclopc calculation suffices to show 
that the spot hypothesis resolves the paradox of granular 
diffusion. Suppose that a spot carries a total free volume, 
Vs, and causes equal displacements, {Axp,Azp), among 
Np particles of volume, Vp. The particle displacement 
can be related to the spot displacement, (Axg, A^;^), by 
an approximate expression of total volume conservation. 




FIG. 7: The Spot Hypothesis II: Each spot displacement 
causes a block of neighboring particles to make a smaller dis- 
placement in the opposite direction, so as to roughly conserve 
the local volume. 



Under this dynamics, the packing of particles is largely 
preserved because a pair of neighboring particles is usu- 
ally affected by the same spot. Occasionally the pair finds 
itself near the edge of a spot, which causes a tiny rela- 
tive displacement. The accumulated effect of such events 
gradually causes the particles to separate as they drain. 
The cage breaking length, however, is much larger than 
the particle size because relative displacements of neigh- 
bors are both small and rare. In fact, the cage-breaking 
length in the model (calculated below) can be compara- 
ble to the system size, as in experiments and simulations 



Ns Vp (Axp, Azp) - -Vs (Ax„ Azp), (10) 

which ignores boundary effects at the edge of the spot. 
Prom this relation, we can easily compute the particle 
diffusion length, 

Var(Axp) w^Var(Ax,) 

o„ = — = — — — 4^ =wb (11) 

^ 2dh\Azp\ 2dhwAzs 

which is smaller than the spot (or volume) diffusion 
length by a factor. 



equal to the ratio of the spot's free volume, Vg, to the 
total affected particle volume, NpVp. The Void Model 
corresponds to the unphysical limit where these volumes 
are both equal to a single particle volume {Np = 1, 14 = 
Vp,w = 1), but generally we expect spots to affect mul- 
tiple particles and carry relatively little excess volume 
{Np ^ l,Vs < IjW <C 1). The nontrivial implication 
of cooperative motion is then that particles diffuse much 
more slowly than volume, bp <C b, which resolves the 
paradox described above. Below we will show that the 
resolution is not just qualitative, but quantitative. 

C. General Formulation of the Spot Model 

The key mathematical concept in the Spot Model is 

that of a diffusing "region of influence" , which causes cor- 
related displacements among all particles within range. 
There are many possible microscopic postulates for this 
local cooperative motion, but the strong tendency to 
preserve nearly jammed packings suggests it should be 
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mainly a block translation as described above, plus per- 
haps a random block rotation, mostly involving nearest 
neighbors. For such rigid-body motions, however, the 
shear at the edge of a spot may be somewhat excessive. 

More generally, we may allow the influence of a spot to 
be governed by a smooth function, u;(rp,rs), which de- 
cays quickly with distance from the spot center, Itj, — r^l. 
For example, in the case of purely translational motion 
in one direction (opposing a spot displacement), the ran- 
dom displacement, AR^'-*, of the ith particle centered at 
rp \ due to the random displacement, AR^-'^ of jth spot 
cented at r^^ would be given by 



AR« 



-7«(r«,rW 



ARl^h arI,^') 



(13) 



We might expect a spot to roughly maintain its "shape" 
while moving through a system of particles, due to the 
local statistical regularity of random packings, so its in- 
fluence function, w;, should depend mainly upon the sep- 
aration vector, — Ts . This may change in regions of 
highly non-uniform flow, but in any case the spot influ- 
ence should typically decay at separations, jr^^ — ri'''|, 
larger than a few particle diameters. Below we demon- 
strate that the mathematical model of Eq. (|13(l is a nat- 
ural starting point for discrete simulations, as well as 
analysis in the continuum limit. 

There are many possible extensions, which should be 
considered in future work. For example, the internal par- 
ticle motion induced by a spot could contain a random 
component, or some other more complicated cooperative 
motion driven by inter-particle forces (see below) . Differ- 
ent spots could also interact with each other by changing 
their drift velocity, diffusion length, size and/or free vol- 
ume, which effectively induces cooperative particle mo- 
tion at large scales. Weak repulsion and attraction might 
both occur, since the volume fraction tends to remain in 
a certain narrow range in a region of flow. Even if spots 
remain unchanged during motion, there could be statis- 
tical distributions of spot free volume, size, and shape, 
which might vary due to compaction and dilation pro- 
cesses driven by random spot annihilation, creation, and 
recombination. Such extensions may be necessary for a 
complete theory of dense granular flow, but in the next 
section we show that the basic model already compares 
well with silo-drainage experiments. 



A. Free Volume and the Peclet Number 

The Spot Model makes a quantitative connection be- 
tween dilation, diffusion, and advection in a dense gran- 
ular flow. As such, it predicts the Peclet number © in 
terms of typical local volume fluctuations. The Peclet 
number is independent of the flow rate, but, unlike with 
the Void Model, the predicted magnitude is consistent 
with experimental data. 

We start with the general observation that the particle 
volume fraction, 4> « 0.6, varies by at most a few percent 
in the region of dense flow, away from the orifice. This 
is certainly a reasonable inference from drainage exper- 
iments, which maintain a nearly uniform density to the 
naked eye js^l . Detailed observations of particles near a 
glass silo wall confirm that the local area fraction varies 
by one to three percent at different fiow rates |lOlj | . Like- 
wise, in granular dynamics simulations, the bulk volume 
fraction in thin horizontal cross sections varies in the 
range, (f> w 0.565 — 0.605, near the orifice, while stay- 
ing within one or two percent of 0.60 in t he region of 
steady dense flow far above the oriflce |l02j |. 

These values are also consistent with studies of random 
packings of hard spheres in zero gravity. Rigidity perco- 
lation and dilatency (expansion under shea r str ess) set 
in at « 0.55 for "random loose packing" |l03l |. which 
provides a rough lower bound on the volume fraction 
in a slow, dense flow. A rough upper bound is set by 
jammed random packing pTj . e.g. in the "maximally 
random jammed state" ( « 0.64 |l3| or (/> « 0.63 [13) 
or the zero-ternperature, zero-stress "jamming point" 
{(f> « 0.63 10,01). Higher volume fractions also pos- 
sible at the expense of some long-range crystalline or- 
der (e.g. as demo nstrated by experiments with horizon- 
tal shaking |104| ). up to the rigorous upper bound of 
(f> = vr/VTS « 0.74 for the FCC lattice 105], but such 
configurations are unlikely to permit quasi-steady fiow. 

If we attribute the maximum local reduction in volume 
fraction, A<j)/(j) « 0.01, to the presence of Ns overlapping 
spots. 



then we have 



1 -I- (l)wNs 



A<j> _ 0.01 



rwNs, 



0.02 



(14) 



(15) 



V. 



EXPERIMENTAL EVIDENCE 



In this section, we consider the simplest version of the 
Spot Model, in which each spot carries the same free vol- 
ume, Vs, and is spherical with a uniform influence, w, out 
to a distance R. We compare this model to experimen- 
tal data, primarily from Ref. |37| . including some new 
results by the same authors. The experiments provide 
compelling tests of the theoretical predictions, because 
the latter were made well in advance llOOll . 



In that case, the Peclet number for horizontal diffusion 
would be 



Pe. 



{\Azp\/At)d _ d 
Var(Axp)/2At ~ 



40 7V, 



(16) 



for b = 1.3d. Since Ng > I in this calculation, we predict 
that a particle falls by on the order of 100 diameters 
before diffusing horizontally by one diameter. 

This simple estimate is consistent with the experimen- 
tal value, PCa: = 321, for particles near a smooth wall in 
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slow drainage from a quasi-two dimensional silo [31]. It 
also suggests that spots occur at fairly high density, with 
as many as Ng — 320/40 = 8 overlapping in a position 
of high dilatency. In a position of low dilatency, where 
A(p/(j> w 0.001, it implies that spots typically do not over- 
lap. Such a limited number of spot overlaps seems like a 
reasonable definition of a "dense" flow, where the simple 
picture in Eq. |7J| could apply. 

(Note that larger Peclet numbers in the range Pe w 
1500 — 3000 have been reported for faster flows in long, 
narrow silos with shear-inducing rough walls 36] . In this 
case, the simple relation in Eq. I|16|) would imply a much 
higher spot density, where positions of 100 overlapping 
spots could be found. At such high spot densities (low 
volume fraction), it seems the model would no longer 
apply due to more independent particle displacements. 
Indeed, since Pe depends weakly on the flow rate, those 
experiments are not in the regime of slow drainage con- 
sidered here. The data is also inconsistent with kinetic 
theories of dilute granular flow [s^, so it seems to cor- 
respond to an intermediate regime of moderately dense 
flow.) 

B. Spatial Velocity Correlations 

The quantitative resolution of the granular-diffusion 
"paradox" provides some support for the Spot Model 
does not directly validate its microscopic hypothesis. 
How can we directly confirm or reject the existence of 
spots? The calculation above suggests that it would be 
impossible to observe the propagation of a single spot, 
but only large numbers of spots, which expectation is 
consistent with x-ray diffraction experiments show ing 
fairly smooth density patterns in draining sand |l06j . 

Rather than trying to observe a single spot, therefore, 
it makes more sense to seek statistical evidence of the 
passage of many spots. A direct signature of the coop- 
erative motion in Fig. [7| is found in the spatial correla- 
tion function of velocity fluctuations (relative to a steady 
mean flow). Two particles are likely to fluctuate in the 
same direction when they are separated by less than a 
spot diameter because most spots engulf them at the 
same time. On the other hand, more distant particles 
are always affected by different spots, which implies in- 
dependent fluctuations. 

The spatial velocity correlation function, C(r), for two 
particles separated by r is easily calculated for uniform, 
spherical spots of radius, R. The two instantaneous par- 
ticle displacements are either identical (perfectly corre- 
lated), if they are caused by the same spot, or indepen- 
dent. Therefore, the correlation function, C{r), is simply 
the scaled intersection volume, a{r; R), of two spheres of 
radius, R, separated by a distance, r: 

°(-«' = '-|s + i^(sr <'^> 

This result appears in a recent study of random point 
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FIG. 8: Direct experimental evidence for spots: Spatial ve- 
locity correlations from the drainage experiments of Ref. Is^ . 
Surface plot (top) and coordinate slices (bottom) of the cor- 
relation coefficient for velocity fluctuations (about the local 
mean flow) at horizontal separation, x, and vertical separa- 
tion, z, in units of particle diameter, d. (Courtesy of Jaehyuk 
Choi.) 



distributions |107| . where a{r;R) is plotted in Fig. 9 
(and earlier in Ref. |l08]). The key point is that C(r) = 
a{r; R) decays to half its maximum value at a separation 
slightly smaller than spot radius, R. Similar curves re- 
sult from other quickly decaying spot influence functions, 
such as exponentials or Gaussians (see below). 

Motivated by this prediction, the experiments in 
Ref. [33 were performed to look for spatial velocity cor- 
relations in a real granular flow. Hundreds of glass beads 
(d = 3mm) near the wall of a quasi-two-dimensional 
draining silo were simultaneously tracked with 1ms res- 
olution in time and d/100 in space in the region of non- 
uniform flow near the orifice. The data clearly reveals 
spatial velocity correlations, shown in Fig. |S1 at the ex- 
pected length scale of several particle diameters. It is 
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interesting that real spots are not quite spherically sym- 
metric. The velocity correlations extend roughly twice as 
far in the horizontal direction as in the vertical direction, 
so these spots are more elliptical in shape. The data also 
implies a smooth decay of the spot influence function, 
w(rp,rs), roughly an exponential decay. In contrast, for 
a uniform spot influence with a cutoff radius, R, the cor- 
relation function, C(r) — a{r,R), would decay roughly 
linearly to zero at r = 2R. 

Leaving these interesting details for future analysis, 
we proceed with our assumption of uniform spot influ- 
ence, as a good first approximation. As expected for 
cooperative motion largely preserving particle cages, the 
spot diameter is roughly five particle diameters, which 
corresponds to tens of affected particles per spot (say 



Np « 20). Using Eq. Hlfci|l to infer spot influence, 



1 



w 



0.0024, 



Pc,6 (321)(1.3) 
the simple relation, Eq. (|12() . then implies 



AT 



0.05 



(18) 



(19) 



In other words, the interstitial free volume carried by a 
spot, although spread across tens of particles, is only a 
few percent of one particle volume. Such small, coherent 
volume fluctuations are completely different from those 
postulated by the Void Model {Vs =Vp, R^ Rp). 



C. Cage Breaking 

Further support for the spot hypothesis comes from 
experimental data on cage breaking j37| , which also pro- 
vides the most direct contradiction of the Void Model. 
Visual observation of particles in draining silos reveals 
that particle cages tend to persist for large distances, 
comparable to the system size. It is not unusual for some 
first neighbors to remain unchanged from the middle of 
the silo all the way to the orifice. A detailed analysis 
of the uniform-flow region far from the orifice yields a 
cage-breaking length, Zc ~ 200(i, extrapolated from small 
falling distances, Az < 5d. As such, the precise value, 
200 d, should not be taken too seriously, but it is clear 
that cage breaking occurs very slowly in regions of low 
shear, at scales comparable to the system size. In con- 
trast, the Void Model inevitably predicts cage-breaking 
at the scale of a single grain, « d. 

Rather than analyzing many-body cage dynamics in 
the Spot Model, it is simpler to consider two tracer par- 
ticles, which are initially first neighbors, in a nearly uni- 
form flow. Without correlations, the horizontal sepa- 
ration, X, would grow with distance dropped, z, like 
X d -\- \J2hiz, where the diffusion length for relative 
motion is twice that of single particles, 62 = 26p, due to 
the additivity of variance. With correlations, the relative 
diffusion length. 



is reduced by a factor, 1 — C(X) k, d/2R (for X k, d). 
Supposing that cage breaking occurs at X « 2d, wc pre- 
dict a cage-breaking length, Zc — d^/262, of roughly. 



d 



R 

2br, 



Pe^R 
2d 



Pe. 



(21) 



where we take the approximation, R w 2d, from the data 
in Fig. |S1 The relation, Zc ~ Pe^d « 300c?, is quite 
consistent with the experimental data. 



D. 



Vertical Diffusion 



We have postulated that spots diffuse upward with a 
horizontal diffusion length equal to the kinetic parame- 
ter, 6^ = b, because the mean velocity roughly follows 
the profile of the Kinematic Model, as shown in Fig. ^ 
By looking at fluctuations, we can also infer the vertical 
diffusion length of spots, which would be difficult to mea- 
sure directly in experiments. (Note that voids exhibit 
zero vertical diffusion, unless one introduces a random 
waiting time between steps.) 

Suppose that the vertical component of the spot dis- 
placement, Azs, is a positive random variable with a fi- 
nite mean, (Az^). We can then define the vertical diffu- 
sion length of spot. 



6 = 



Ya,r{Azs) 
2dh{Az,) 



(22) 



by analogy with Eq. ^ for the horizontal diffusion length 
(where Az^ should be replaced by (Az^)). In a similar 
way, we define the vertical diffusion length of particles, 
bl (and b^ = bp). 

Assuming that the spot inffucnce does not distinguish 
between different velocity components, as in Eq. H10() . we 
have the identity, 



bj,^b_^ 
hW b^ 



(23) 



This allows us to infer the vertical diffusion length of 
spots from that of particles. 



b^l^'l^ 



w Pe^w (150)(0.0024) 



2.8 d 



(24) 



b2{X)^2bp[l-C{X)], 



(20) 



where we take Fez = 150 from experiments [23. This cal- 
culation assumes that the dominant contribution to the 
particle vertical diffusivity is the influence of the spot 
vertical diffusivity, and not the random arrival of dif- 
ferent spots (consistent with the mean-field analysis of 

section rVIl|) . Since bl/bj- = 2.8/1.3 w 2, the spot dis- 
placements in these experiments fiuctuate twice as much 
in the vertical direction as in the horizontal direction. 
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avoid an accident by a flux-density relation, q{p), with a 
maximum at pc- The conservation law for the car density, 



FIG. 9: Sketch of nonlinear waves of fre e vo lume (gray re- 
eions) observed in funnels of draining sand For the same 

orifice width, the waves propagate upward for large opening 
angles (a) and sometimes downward for small opening angles 
(b). 



E. Spot Interactions and Density Waves 

Until now, we have assumed that spots undergo inde- 
pendent random walks. This suffices to describe many 
aspects of silo drainage, so we conclude that spot inter- 
actions are relatively weak in such flows. Nevertheless, 
there must be short-range interactions between spots to 
keep the local volume fraction in the typical range for 
dense, gravity-driven flows (0.56 < </> < 0.61). Physically, 
spot repulsion should derive from the tendency for parti- 
cles to quickly collapse from all directions toward a region 
of overly low volume fraction {(p < 0.56). Because each 
spot carries very little volume, the repulsion should only 
become strong when many spots overlap (e.g. Ng > 10). 
Similarly, spot attraction follows from the tendency of 
free volume to avoid nearly jammed regions (</> > 0.61) 
if it can be accommodated in an already flowing region. 
This attraction should strengthen as spots become more 
dilute {Ns < 1), e.g. sharpening the boundary between 
flowing and stagnant regions in Fig. 

Is there any evidence for spot interactions in exper- 
iments? An indirect sign could be the density waves 
observed by Baxter et al. in three-dimensional drain- 
ing funnels of rough sand using x-ray diffraction |l06l |. 
The authors gave a heuristic explanation of such waves 
in terms of the formation and collapse of arches across the 
central region but did not attempt a quantitative theory. 
As sketched in Fig. O the observed density patterns were 
mostly consistent with upward moving waves of reduced 
density, but downward moving waves were also inferred 
in some narrow funnels. Nonlinear wave equations for di- 
lute flows 1 l)!f have been derived from the kinetic theory 
of inelastic gases , but there does not appear to be a 
microscopic theory for density waves in dense flows. 

Subsequent authors jTlS IHl [ill have suggested 
an analo gy with nonlinear kinematic waves in traffic 
flow |ll3l | , and some have even claimed that granula r flow 
and traffic flow are in the same "universality class" |ll4j . 
The classical, continuum theory of traffic waves describes 
the effective short-range repulsion between cars trying to 



dp dq dp dp 



(25) 



then involves the characteristic velocity, c(p) = q'{p), 
which changes sign at pc, implying that density waves 
propagate forward in light traffic (p < pc) and backward 
in heavy traffic (p > Pc)- Whenever density variations 
are generated (by noise, boundary conditions, etc.), they 
combine into coherent structures terminated on one side 
by a shock and on the other by an expansion fan. 

Granular shock waves in two-dimensional inclined fun- 
nels and pipe s have bee n studied extensively in sim- 
ulation s [lUl Ill2l Ill5l | and particle-tracking experi- 
ments |62l l63l| . A detailed study of shock statistics 
has revealed differences with simple traffic-flow models, 
although somewhat similar average q{p) curves have been 
measured. Granular dynamics simulations also show that 
friction among the particles and walls initiates the fluc- 
tuations leading to shock formation IIOJ. 

The Spot Model provides a simple microscopic expla- 
nation for the experimental observations. Given our ex- 
pectation of weak spot repulsion, it is natural to think 
of spots as "cars" moving upward in a dense flow. For a 
given orifice width, a funnel with a large opening angle, 
as in Fig. I^a), should have a fairly low spot density as 
a result of lateral diffusion. In that case (p^ < pc), non- 
linear spot waves initiated by ffuctuations at the orifice 
propagate upward at nearly the spot drift velocity, slowly 
organizing into coherent structures of reduced density. 
The backside of the wave forms shock, below which the 
particle density suddenly increases. Particles gradually 
accelerate into the looser spot-rich region and suddenly 
slow down upon entering the more compact region below. 

As the funnel angle is reduced (for the same orifice 
width), spots are increasingly reffected by the bound- 
aries into the central region, where their density goes up. 
The crowding of spots could eventually cause a reversal 
of density waves, analogous to the familiar backward- 
moving waves of heavy traffic. In that case (ps > pc), 
a region of reduced spot density would form a shock on 
its back (upper) side, which the particles would see as 
downward moving compaction wave. On the other hand, 
a region of increased spot density would form a shock 
on its front (lower) side, corresponding to a downward 
moving rupture zone. 

It would be interesting to revisit experimental density 
patterns in light of the Spot Model. Density waves in 
wide funnels (ps < Pc) should move as weakly interact- 
ing "tracer spots", thus reveahng the spatial profile of 
spot drift and diffusivity. For example, the wave patterns 
sketched in Fig. Ela) suggest that the spot drift velocity 
and diffusion length are smaller near the stagnant region 
than in the faster-flowing central region. 
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F. Crossover from Super to Normal Diffusion 

We have already mentioned the surprising result of 
particle tracking in draining silos that all fluctuations 
depend on the distance dropped, not time, over a wide 
range of flow rates Wi\ ■ For the mean-squared displace- 
ment relative to the mean flow, there is a universal 
crossover from super-diffusion to normal diffusion after 
a particle falls by roughly one diameter. The observed 
super-diffusive regime, O.OOSd < |Aa:| < 0.5d, is not bal- 
listic (|Aa;| oc t (X z) but instead exhibits a non-trivial 
scaling exponent. Ax oc t"-^^ c>c z"-^^. 

A small-scale ballistic regime has also been inferred by 
diffusing wave spectroscopy for glass beads {d = 95/xm, 
V = 0.32 c m/s) flowing down a vertical pipe toward a 
wire mesh jllfif . The data does not seem of indicative 
of thermal collisions, however, since the "collision dis- 
tance" , Ic — "28 nm « 10~^(i, is surely at the scale of sur- 
face roughness. Moreover, the "collision time" , Tc = 9 /iS 
« 10~'^d/v, implies a typical velocity, Ic/tc = 0.31 cm/s, 
equal to the mean flow speed, which seems hard to at- 
tribute to independent thermal fluctuations. In any case, 
no transition to diffusion is observed at this scale, so these 
nano-coUisions are not of the randomizing type postu- 
lated by kinetic theory and have little bearing on our 
discussion of structural rearrangements. 

Returning to the particle-tracking data, we have seen 
that the Spot Model describes the diffusive regime with 
quantitative accuracy, but the super-diffusive regime is 
not so easily explained. Spots are not constrained to 
move on a lattice, so the size of the spot displacements is 
arbitrary, as long as the horizontal spot diffusion length 
equals the kinetic parameter, = h ~ d. For example, 
each spot could follow a persistent random walk in which 
successive, small random displacements (Az^ <C d) have 
a correlation coefficient, 7. Following Taylor |91j|,|93|, the 
kinetic parameter would then be 



where 60 is the spot diffusion length in the absence of 
auto-correlations, given by Eq. 1^. For < 7 < 1, the 
persistent random walk of the spot has ballistic scaling 
up to a crossover distance, —Azg/ log 7, which should be 
of order d, since b ~ d. Since the spot diameter is larger 
than d, a single particle is affected many times by the 
same spot, so the particle would inherit a persistent ran- 
dom walk with ballistic motion up to a crossover distance 
of roughly w d <^ d. 

This scaling and crossover distance do not agree with 
experiment, but various modiflcations could perhaps re- 
solve the discrepancy. It is straightforward to achieve 
sub-ballistic scaling with longer-range, power-law auto- 
correlations between spot steps, although some cutoff 
would be needed to preserve the normal diffusive regime. 
More likely, however, the super-diffusive regime is asso- 
ciated with more complicated particle motion caused by 



spots. The simple picture in Fig. |7| ignores packing con- 
straints and frictional contacts, so we should expect ad- 
ditional microscopic rearrangements to occur, superim- 
posed on the mean spot-induced motion. This funda- 
mental issue becomes more clear in simulations of the 
model. 



VI. SPOT-BASED SIMULATIONS 

A. Statistical Dynamics of Dense Flow 

The Spot Model provides a very simple and efflcient 
Monte Carlo algorithm for dense granular drainage. The 
simulation begins with a given distribution of tracer par- 
ticles, either on a lattice or in a random packing. The 
particles then move passively in response to spots under- 
going directed random walks upward. Spots are injected 
randomly in time and space along the orifice to artificially 
set the flow rate. (Our theory makes no prediction about 
how spots are created, which is related to the poorly un- 
derstood connection between flow rate and oriflce width.) 

For illustration purposes, we consider the simplest ver- 
sion of the model in two dimensions without boundaries. 
Let us assume that spots undergo independent, directed 
random walks upward on a lattice (without interactions). 
The vertical spot lattice parameter is set on the order of a 
particle diameter, d, and the horizontal lattice parameter 
is constrained by kinetic parameter, 6, via Eq. I^. Alter- 
natively, for smoother spot motion using a persistent ran- 
dom walk, the vertical lattice parameter may be reduced 
for a desired level of auto-correlation, say 7 = 0.5, where 
b is now set by Eq. H26|l . We further assume Eq. H10() for 
the spot influence, where each particle within a distance, 
i?, from the spot center (after a displacement), receives 
the same anti-parallel displacement, smaller by a factor, 
w. 

Spot simulations of the experiments in Fig. and [51 
are shown in Fig. ^| and ^] respectively. The spot 
parameters determined from independent experimental 
measurements, as described above {b = 1.3 d, R = 2.5 d, 
w — 0.0024), without any ad hoc fitting. The result is a 
dramatic improvement over the simulations in Figs.0]and 
ISJc) of the same situations using the Void Model, which 
produces far more mixing. Remarkably, the mean flow 
proflle is nearly the same in all cases, given by Eq. 
even though the particle dynamics is radically different. 

It is visually apparent that the Spot Model resolves 
the "paradox of granular diffusion" . Detailed analysis 
(to appear elsewhere) shows that the experimental dy- 
namics of two-color interfaces is reproduced rather well 
by simulations, as in Fig. 1111 The interface remains sharp 
while moving with the mean flow and slowly developing 
fluctuations at the scale of the spot size. 
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FIG. 10: A spot simulation of the experiment in Fig. to 
be compared with the void simulation in Fig.|l] (Courtesy of 
Jaehyuk Choi.) 
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FIG. 11: A spot-model simulation (right) of the experiment 
in Fig. El (left). 



B. Density Fluctuations and Packing Constraints 

The reader may have noticed something strange about 
the simulations in Figs. llUl and ITTI — the nearly uniform 
density of particles is not preserved. In regions of shear, 
the particles form dense overlapping bands separated by 
thin vacuum bands transverse to the shear direction. Al- 
though mesoscopic shear bands do form in disordered 
materials, such as metallic glasses (see section Hxj) . these 
bands are at the particle scale, and hence are not physi- 
cal. At the lower edges of the flow region, where the shear 
rate is highest, the fluctuations intensify, while directly 
above the orifice the density increases. 

We conclude that the simplest version of the Spot 
Model does not completely describe the joint probability 
distribution of all particles, even though it predicts the 
marginal distribution of each tracer particle quite well. 
As expected, what is missing is a strict enforcement of 
packing constraints. The block-like cooperative motion 
in spots causes density fluctuations to grow very slowly, 
but unphysical structures eventually arise. 

Stabilizing the density with a simple analytical mod- 
ification is likely to be difficult. As discussed above, 
density-preserving extra fiuctuations appear to be re- 
lated to the super-diffusive regime seen in experiments, 
which has a nontrivial, sub-ballistic scaling. Moreover, 
we should not expect to be able to generate a multitude 
of random close packings by a simple statistical model. 
Even without flow, the sampling of random close pack- 
ings is a challenging computational problem in its own 
right, which seems to req uire brute- force molecular dy- 
namics [Ulillillillillli 

On the other hand, unphysical density fiuctuations also 
arise in the standard kinetic theory of liquids, when parti- 
cles undergo independent random walks. In dilute gases, 
this assumption correctly predicts Poisson-distributed 
density fluctuations, where the v arian ce in the particle 
number scales with the volume |117| . but in a liquid, 
just as in a granular material, this picture is incorrect. 
Packing constraints during flow somehow preserve much 
smaller, "hyper-uniform" density fluctuations, where the 
variance scales with the surface area, as in the case of a 
regular lattice 107] . It is fair to say that a collective sta- 
tistical theory of liquids, molecular or granular, has not 
yet been developed, so we should not be easily deterred. 



C. Multiscale Algorithms for Dense Flows 

We have seen that the collective particle dynamics in 
the Spot Model is quite realistic, aside from slowly grow- 
ing density fluctuations, so it may be that only occasional 
small perturbations could suffice to preserve the density. 
For example, after some (or all) spot-induced cooperative 
displacements, the affected particles could be relaxed to 
equilibrium with simple spring-like forces, perhaps in- 
cluding the effect of gravity. The center of mass of the 
block of displaced particles should be constrained during 
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FIG. 12: The Spot Hypothesis III: Real spots also cause small 
internal rearrangements driven by packing constraints (rela- 
tive to the average motion shown in Fig. |7J . In a multiscale 
simulation algorithm, this could be done by a constrained re- 
laxation following the mean spot-induced displacements. 



relaxation to prevent slipping back into the initial posi- 
tions. A preliminary implementation of this idea using 
a soft-core repulsion between grains (without gravity or 
friction) prod uces rather realistic bulk flows with a stable 
density |lTal |. 

Regardless of the detailed implementation, fixing the 
density will result in some additional cooperative fluctu- 
ations relative to the mean spot-induced displacement, 
as shown in Fig. 1121 This is surely a more realistic pic- 
ture of a dense flow, which incorporates the microscopic 
intcr-particle forces into a new kind of multiscale sim- 
ulation algorithm. In spite of the potentially expensive 
relaxation step, such an algorithm is much more efficient 
than brute-force molecular dynamics tracking all parti- 
cle contacts. It takes advantage of spots to accelerate the 
dynamics, while keeping the system in a nearly jammed 
state. 

This would be an interesting direction to pursue, not 
only for simulations, but also to gain insight into the ba- 
sic physics of dense flows. Unlike the mean spot-induced 
motion, the extra fluctuations to preserve the density are 
correlated over the passage of multiple spots. The result 
could be a crossover from super to normal diffusion af- 
ter a particle interacts with many spots, perhaps once it 
falls by its own diameter. Understanding such small-scale 
fluctuations is a fundamental open question. 

VII. THE CONTINUUM LIMIT 

Macroscopic approximations offer many advantages 
over microscopic theories, whenever a clear connection 
can be made. Both analytical and numerical solutions 
are usually much easier starting from a continuum ap- 



proximation than an underlying discrete model. More- 
over, various poorly understood microscopic details are 
swept into parameters in the continuum equations, which 
(one hopes) capture the essential physics in a robust way. 
A good example is the simultaneous description of Brow- 
nian motion, flnancial time series, and heat transfer by 
the diffusion equation. Perhaps even more remarkable is 
continuum hydrodynamics, which hides our ignorance of 
collective statistical dynamics in liquids. In this section, 
we derive continuum equations for diffusion and flow in 
the Spot Model based on the simple picture of Fig.[7| thus 
ignoring the complicated density-preserving relaxation in 
Fig.im 

A. A Nonlocal Stochastic Differential Equation 

We begin by partitioning space as shown Fig.^l where 
the nth volume element, AV^i"', centered at ri"' contains 
a random number, AA^i"\ of spots at time t (typically 
one or zero). In a time interval. At, suppose that the jth 
spot in the nth volume element makes a random displace- 
ment, AR^-'' (ri"-*) (which could be zero). According to 
Eq. l(T^ . the total displacement, ARp, of a particle at Tp 
in time At is then given by a sum of all the displacements 
induced by nearby spots, 

ARp(rp) = -^ z«(rp,r(")+AR(^-)(r:)) AR5")(r:) 

n j — 1 

(27) 

Note that the spatio-temporal distribution of spots, 
AA^i"', is another source of randomness, in addition to 
the individual spot displacements, AR^"-*, so that each 
particle displacement is given by a random sum of ran- 
dom variables. 

In the continuum limit, we arrive at a non-local, non- 
linear stochastic differential equation (SDE): 

rfRp(rp,t) = - J dNs{rs,t)w{rp,rs+d'Rs{rs,t))dIis{rs,t), 

(28) 

where the stochastic integral is defined by the limit of the 
random Riemann sum in Eq. I|27|) . This equation differs 
from standard nonlinear SDEs in two basic ways: (i) The 
tracer trajectory, 

rp(t)= / rfRp(rp,T) (29) 

Jt=0 

is passively driven by a stochastic distribution of mov- 
ing infiuences (spots), dNs{rs,t), which evolves in time 
and space, rather than by some internal source of in- 
dependent fluctuations, and (ii) the stochastic differen- 
tial, c?Rp(rp, t), is given by a non-local integral over other 
stochastic differentials, c?Rs(rs,t), associated with these 
moving influences, which lie at flnite distances away from 
the particle. 
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and difFusivity tensor, 



(r,i) 



2dt 



(32) 



dRJr;) 




dRJrJ 



dN= 



dN= 1 



FIG. 13: Sketch of a particle interacting with a coUection 
of passing spots, showing some of the quantities involved in 
defining the nonlocal SDE, Eq. ij^ . 



B. Mean- field Fokker-Planck Equation 

Equation H28|l seems to be a new type of SDE, for 
which a mathematical theory needs to be developed. 
Here, we begin deriving analytical results based on the as- 
sumption that all the stochastic differentials for particle 
displacements in Eq. H29() are independent. This ignores 
the fact that the spot distribution at one time, dNs{rs, t), 
depends explicitly on the distribution and displacements 
at the previous time, t — dt. Even for a stationary spot 
distribution, it also ignores the fact that each spot, due 
to its extended influence, generally affects the same par- 
ticle more than once, so that any autocorrelations in spot 
motion are partly transferred to the particles. 

Assuming independent displacements, the propagator, 
Pp{r, t\ro, to), which gives the probability density of find- 
ing the particle at r at time t after being at Tq at time to, 
satis fies a Fokker-Planck (or forward Kolmogorov) equa- 
tion [TT^ . 



dPp 
dt 



V • (upPp) = VV 



(30) 



with drift velocity. 



Up(r,t) 



dt 



j.^^ (ARp(r,t)) 

At^O At ' 



(31) 



(Here VV : A denotes J^c^Ep w§i;-) The Fokker- 
Planck coefficients may be calculated by taking the 
appropriate expectations using Eq. H27|) in the limits 

Ay/"' and then At 0, which is straightforward 
if we assume (as above) that spots do not interact. In 
this approximation, the spot displacements, AR^^'(ri"^), 
and the local numbers of spots, ATVi"-*, are independent 
random variables within each infinitesimal time interval. 
They are also independent of the same variables at earlier 
times, as assumed above. 

In order to calculate the drift velocity, we need only 
the mean spot density, ps(rs,t), defined by (AA^i"') — 
Ps{rs,t) AVs^"\ The result, 

Up(rp,t) = - J dVsw{rp,rs) [ps(rs, t)us(rs, t) 

~2-Ds{rs,t)-Vps{rs,t)] (33) 

exhibits two sources of drift. The first term in the inte- 
grand simply opposes the spot drift velocity. 



Us(r,i) 



(rfRs(r,t)) 
dt 



(34) 



as in the heuristic equation (|10l) . The second term, which 
depends on the spot diffusion tensor. 



^ 2dt ' 



(35) 



is a "noi se- induced drift" (typical of nonlinear 
SDEs |119| |') causing particles to climb gradients in 
the spot density. Both effects are averaged over 
nearby regions, weighted by the spot influence function, 
w{rp,rs). 

In order to calculate the diffusivity tensor, we also need 
information about fluctuations in the spot density. Con- 
tinuing with our "mean-field approximation" , it is natu- 
ral to assume independent spot fluctuations, 

(AAi'^^AAi")) = <5™^„((AAi"))2) = O ((AF^"))"^) , 

(36) 

where v ~ I for a Po isson process and i> < 1 for a hyper- 
uniform process |107| . It turns out that such fluctuations 
do not contribute to the diffusion tensor (in more than 
one dimension), and the result is 

Dp(rp,t)== j dVsw{Yp,Ysfps{Ys,t)-D,{vs,t). (37) 

Note that the influence function, w, appears squared in 
Eq. (|37|l and linearly in Eq. (|33|l . This causes the Peclet 
number for tracer particles to be of order 1/w smaller 
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than that of spots (or free volume), as was derived heuris- 
tically above in Eq. 1)11(1 . 

An interesting general observation about Equation 
(|30|l with coefficients and l|37|l is that rescaling the 
spot density is equivalent to rescaling time. In other 
words, all aspects of the stochastic motion of a tracer 
particle, including its Peclet number, are determined by 
geometry, independent of the total flow rate (or average 
spot density). As explained above, this experimental ob- 
servation is in qualitative disagreement with kinetic the- 
ory and quantitative disagreement with the Void Model. 
On the other hand, the fact that fluctuations are pro- 
portional to the spot density suggests that a meaningful 
definition of "temperature" for dense granular flows (if 
one exists) might simply depend on the free volume of 
random packings. This possibility is discussed below in 
section Hxl 

Higher-order terms a Kramers-Moyall expansion gen- 
eralizing Eq. (|30|l for finite independent displacements, 
which depend on fluctuations in the spot density, are 
straightforward to calculate, but beyond the scope of this 
paper. Such terms are typically ignored in stochastic 
analyses because, in spite of improving the overall ap- 
proximation, they tend to produce small negative proba- 
bilities in the tails of distributions jll9l | . For dense granu- 
lar flows described by Eq. (|28l) . however, they may be im- 
portant since displacements can be relatively large com- 
pared to spatial variations in probability density. This 
is the case in highly nonuniform flows, such as near the 
orifice of a silo, where mean velocities may vary on the 
scale of a few particle diameters. 



C. Spatial Velocity Correlation Tensor 

For any stochastic process representing the motion of 
a single particle, it is well-known that transport coeffi- 
cients can be expressed in terms of tem pora l correlation 
functions via the Green-Kubo relations |ll9f . For exam- 
ple, the diffusivity tensor in a uniform flow is given by 
the time integral of the velocity auto-correlation tensor, 



p 



(38) 



where Vp{t) — {Up} — dRp/dt is the stochastic velocity 
of a particle. (A similar relation holds for spots.) 

In the Spot Model, nearby particles move coopera- 
tively, so the transport properties of the collective system 
also depend on the two-point spatial velocity correlation 
tensor, 



^{U»{v,Y){U^{v^) 



(39) 



which is normalized so that Cp^{r, r) = 1. We emphasize 
that the expectation above is conditional on finding two 
particles at ri and r2 at a given moment in time and 



includes averaging over all possible spot distributions and 
displacements. Substituting the SDE into Eq. JSHJ 
yields 

^cpt , _ I dVs Psjrs) Dfjvs) w{v,,Vs) wjv^.Vs) 
Op ^ri,r2j— ; 



(40) 

where we have assumed again that spots do not interact 
(independent displacements). 

Equation (|40|l is a new kind of integral relation for 
cooperative diffusion, which relates the spatial velocity 
correlation tensor to the spot (or free volume) diffusiv- 
ity tensor via integrals of the spot influence function, 
w{vp,Ys)- If the statistical dynamics of spots is homoge- 
neous (in particular, if is constant), then the relation 
simplifies: 



C;'^(ri,r2) 



/ dVs Ps{Ys)w{yi,Ys)w{y2,Ys) 



dVs Ps{ys) u;(ri, r,)2 / dV^ p,(rO w{Y2,Y',y 

(41) 

if Df^ 7^ (and otherwise). If the statistical dynamics 
of particles is also homogeneous, as in a uniform flow 
(ps — constant), then it simplifies even further: 



„/3/„x ^ JdVsW{Y-Y,)w{~Y,) 



(42) 



where we have assumed that the spot influence function, 
and thus the correlation tensor, are translationally invari- 
ant (r = ri — r2). These predictions may be used to infer 
the spot influence function experimental measurements, 
as in Eq. H17I) above for uniform spots with sharp cutoff. 
For the case of a symmetric Gaussian spot of width a in 
each direction. 



u;(r) 



Vs 



(43) 



(27rcr2)3/2 

an even simpler formula results from the overlap integral 



Cf{Y)=e 



(44) 



D. 



Relative Diffusion of Two Tracers 



The spatial velocity correlation function affects many- 
body transport properties. For example, the relative dis- 
placement of two tracer particles, r = ri — r2, has an 
associated diffusivity tensor given by. 



Z?"^(ri,r2) = i?^'^(ri)-f i?^'^(r2) 



— n"/3 



a/3/ 



(45) 



2Cf{Y,,Y,)jD^-{Y,)D^pP{Y,) 



In a uniform flow, the diagonal components take the sim- 
ple form 



D°"^{y) =2D'^'^ {l-C^'^iY)) 



(46) 
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which was used above to estimate the cage-breaking 
length. A more detailed calculation of the relative propa- 
gator, P{r,t\ro,to), neglecting temporal correlations (as 
above) would start from the associated Fokker-Planck 
equation, 



dPjr, t) 
dt 



VV : (D(r)P(r,i)) 



(47) 



with a delta-function initial condition. (In a non-uniform 
flow, one must also account for spurious drift and motion 
of the the center of mass.) 

Our simplified analysis does not enforce packing con- 
straints, so it allows for two particles to be separated 
by less than one diameter. A hard-sphere repulsion may 
be approximated by a reflecting boundary condition at 
|r| = d when solving equations such as H47|l . but there 
does not seem to be any simple way to enforce inter- 
particle forces exactly. However, this may be precisely 
why the model remains mathematically tractable and in- 
tuitive, while capturing the essential physics of diffusion 
in dense flows. 



VIII. APPLICATION TO SILO DRAINAGE 
A. Statistical Dynamics of Spots 

The analysis in the previous section makes no assump- 
tions about spots, other than the existence of well-defined 
local mean density, mean velocity, and diffusion tensor, 
which may depend on time and space. As such, the re- 
sults may have relevance for a variety of dense disordered 
systems exhibiting cooperative diffusion (see below). In 
this section, we apply the general theory to the specific 
case of granular drainage, in which spots diffuse upward 
from a silo orifice. 

For simplicity, let us assume that each spot undergoes 
mathematical Brownian motion with a vertical drift ve- 
locity. Us = VgZ, and a diagonal diffusion tensor. 




(48) 



which allows for a different diffusivity in the horizon- 
tal (_L) and vertical (||) directions due to the symmetry- 
breaking effect of gravity. In that case, the propagator 
for a single "spot tracer" , Ps(x, z, t|xo, zq, i), satisfies an- 
other Fokker-Planck equation. 
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(49) 



The coefficients may depend on space (e.g. larger velocity 
above the orifice than near the stagnant region), as sug- 
gested by the e xperi mental density-wave measurements 
discussed above |l06j . 

The geometrical spot propagator, ^'^(x, jz, Xq, zo), is 
the conditional probability of finding a spot at horizontal 



position X once it has risen to a height z from an initial 
position (xq, zq). For constant Vg and Dg, the geometrical 
propagator satisfies the diffusion equation. 



dVs 
dz 



(50) 



where b = D-^Vg is the kinematic parameter. If spots 
move independently, this equation is also satisfied by the 
steady-state mean spot density, ps(x, z), analogous to 
Eq. 0) for the void density in the Void Model. However, 
the mean particle velocity in the Spot Model, Eq. H33|l. 
is somewhat different, as it involves nonlocal effects (see 
below). 

The time-dependent mean density of spots, /3s(x, z, t), 
depends on the mean spot injection rate, Q{xo,zo,t) 
(number/areaxtime), which may vary in time and space 
due to complicated effects such as arching and jamming 
near the orifice. It is natural to assume that spots are 
injected at random points along the orifice (where they 
fit) according to a space-time Poisson process with mean 
rate, Q. In that case, if spots do not interact, the spatial 
distribution of spots within the silo at time t is also a 
Poisson process with mean density, 

Ps{x,z,t) = dxf) dzf) / dto 

J J Jta<t 

Q(xo, zo, t) Ps{x., z, i|xo, zo, h)- (51) 

For a point-source of spots (i.e. an orifice roughly one 
spot wide) at the origin with fiow rate, Qo{t) (num- 
ber/time), this reduces to 



/9,(x,z,i)= / dioQo(io)n(x,z,t|0,0,to), (52) 

Jto<t 

where P, is the usual Gaussian propagator for Eq. (|49|l in 
the case of constant Us and D^. In reality, spots should 
weakly interact, but the success of the Kinematic Model 
suggests that spots diffuse independently as a good first 
approximation. 



B. Statistical Dynamics of Particles 

Integral formulae for the drift velocity and diffusivity 
tensor of a tracer particle may be obtained by substitut- 
ing the spot density which solves Eq. (|51|l into the general 
expressions (|33|) and (|37|l . respectively. For example, if 

spots only diffuse horizontally (Z)| = 0), then the mean 
downward velocity of particles is given by 



Vp{r,t) ^ / dVsw{rp,rs) Ps{rs,t)vs{rs,t) 



(53) 



The main difference with the analogous prediction of the 
Void Model, Eq. is the convolution with the spot 

influence function. These sorts of integrals are most im- 
portant in regions of high shear, where the spot density 
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and/or spot dynamics varies on scales comparable the 
spot size (several particles wide). 

For simplicity, let us consider a bulk region where the 
spot density varies on scales much larger than the spot 
size. In this limit, the integrals over the spot influence 
function reduce to the following "interaction volumes" : 



(54) 



for k = 1,2. (Note that Vi = Vg above.) The equation 
for tracer-particle dynamics H3U|) then takes the form, 



dPp 
dt 



d_ 

dz 



VsP. 



2D 



\dps 



dz 



ViP„ 



(dIpsV^Pp) + Vi {D^PsV^Pp) . (55) 



Again, it is clear that rescaling the spot density is equiv- 
alent to rescaling time. 

When the spot dynamics is homogeneous (i.e. Us and 
T)s are constants). Equation H55I) simplifies further: 



1 dPp 
dt 



dz 

2b^V ■ (PpVps) 



^2 

Kq;2 ) (p^Pp) 



(56) 



2b\\^(p^ 
dz V' dz 



where b 



b = D^/vs and 6" = d\/vs are the spot 
diffusion lengths and b^ = bpV2/Vi and b^, — 6IIV2/V1 
are the particle diffusion lengths. In this approximation, 
the latter are given by the simple formula, 



p 

b^ 



4 _ / dVs w{v, Vs 



bW JdVsw{r,rs) 



(57) 



which generalizes Eqs. Hll|) and (|23l) derived above for 
a uniform spot with a sharp cutoff. In the case of a 
Gaussian spot, Eq. H43|l . we obtain a simple relation. 



b^ foil 87r3/2cr3 



(58) 



between the particle and spot diffusion lengths and the 
spot size and free volume. 

The physical meaning of the diffusion lengths becomes 
more clear in the limit of uniform flow, ps = constant. In 
terms of the position in a frame moving with the mean 
flow, C = Vpt — z, where Vp = VsVsPs, we arrive at a 
simple diffusion equation, 



dPp 
dC 



bi 
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Pdz^ 



P. 



(59) 



where C acts like time. Therefore, bp and bp are the 
variances of the displacements in the perpendicular and 
parallel directions per twice the mean distance dropped. 

Solutions to the new equations above and comparisons 
with experimental data are left for future work. 



IX. OUTLOOK FOR GRANULAR MATERIALS 
A. Silo Drainage as a Physics Problem 

Gravity-driven silo drainage is an ideal system for fun- 
damental studies of diffusion and mixing. The simplic- 
ity of the setup, free of externally applied shear, rota- 
tion, or vibration, may make it the easiest to under- 
stand in microscopic detail. Future work should focus on 
detailed comparisons between the present theory (both 
discrete models and continuum equations), experiments, 
and granular dynamics simulations. It is natural to start 
with dry, cohesionless, monodisperse grains, as in this 
work, and then move on to more complicated granular 
materials. Surely many more physical principles remain 
to be discovered. An important goal should be connect 
with an appropriate continuum model, perhaps among 
those cited in the introduction. Few of these models have 
been carefully tested in silo-drainage experiments, and 
none has been derived systematically from a microscopic 
theory including diffusion. 

One such model, due to Aranson et al. [s^ [s^ pos- 
tulates the co-existence of "liquid-like" and "solid-like" 
micro phases in the granular flow, where the local liquid 
fraction is an order parameter following ad hoc Landau- 
type dynamics. (This is a substantial generalization of 
models of surface avalanches, which postulate co-existing 
static and rolling phases 120, 121], and bear s som e sim- 
ilarity to models of the glass transition j ^sl Il22| .') The 
length scale for order-parameter variations is assumed to 
be of order the particle diameter, even though this might 
seem too small for a continuum theory of co-cxisting 
phases. 

In spite of its successes in describing experimental 
flows, the order parameter lacks any clear microscopic 
basis. The present work, however, suggests that it may 
be related to the spot density, at least in the case of 
silo drainage (which has not yet been analyzed). If a 
successful connection could be made with statistical dy- 
namics, it would be reminiscent of the microscopic theory 
of superconductivity, which came decades after the phe- 
nomenological macroscopic description of Ginzburg and 
Landau ^58] . Such connections may also exist with other 
empirical continuum theories of granular and glassy sys- 
tems, which invoke the "free volume per atom" as a sort 
of order parameter controlling material response. 



B. Some Open Questions about Spots 

There are many questions about silo drainage raised 
by this work, once one accepts the existence of spots. 
How do spots move? Does their drift velocity, diffusivity, 
and/or shape depend on particle properties, position in 
the silo, or the presence of other spots? How do spots 
interact with container walls or free surfaces? How do 
particles move within spots? What are the "extra" fluc- 
tuations associated with small-scale super-diffusion and 
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the strict enforcement of packing constraints? Does the 
spot influence depend on the local velocity, shear rate, 
or stress, the friction and elasticity of the particles, etc.? 
How do spots interact? What kind of attraction and re- 
pulsion between spots is responsible for stabilizing the 
volume fraction in the range of flowing random pack- 
ings? Do spots subdivide and recombine? How are spots 
created and destroyed? How do spots behave at an orifice 
or free surface? How do dilation and compaction create 
and destroy spots in the bulk? 

Even more basic questions have to do with the founda- 
tions of the theory. Why do spots even exist? In a general 
sense, cooperative motion arises from strong dissipation 
in the presence of packing constraints (as in gravitational 
inelastic collapse [lO, Il23j ) , but why should the motion 
be driven by coherent spots of free volume as described 
here? This motivates the question of whether spot prop- 
erties can be predicted from first principles. If so, what 
are the key ingredients of a complete ab initio theory? It 
seems that spots are mainly a consequence of geometrical 
constraints in flowing, dense random packings, but this 
connection should be made more precise. The statistical 
geometry of random packings may be an appropriate 
starting point. 



C. Other Gravity Driven Flows 

Since there is strong evidence for spots in silo drainage, 
it is natural to assume they also exist in other gravity- 
driven flows. In the case of flow down an inclined 
chute |124| . cooperative motion has also been invoked, 
albeit without a precise microscopic description. Erta§ 
and Halsey have explained the Bagnold scaling of the 
mean velocity in terms of "granular eddies" : coher- 
ent structures, several particles in diameter, which ro- 
tate ri gidly in the same sense (as if rolling down the 
plane) |l23l |. However, no evidence for granular eddies 
has yet bee n fou nd in experiments or granular dynamics 
simulations |102| . 

The basic argument for granular eddies Il23( could also 
be made for spots: Strong dissipation ("gravitational col- 
lapse" causes particles to move cooperatively with 
their nearest neighbors. Perhaps some composite picture 
is appropriate for chute flows, although internal rotation 
with the shear may not be necessary. With simple spots, 
shear can occur, as in the flow region in a draining silo, 
when there is a gradient in the spot density transverse to 
the mean spot drift. With a theory of how spots interact 
with the free surface, it may be possible to explain the 
data without including internal same-sense rotations (or 
eddies) in the Spot Model. 

In the context of similar experiments in an inclined 
pipe [s^ , Pouliquen et al. have proposed a totally differ- 
ent mechanism for dense shear flow j53] . They postulate 
that flow occurs by a nonlocal fluctuation- activated pro- 
cess: A shear event in one location initiates stress fluc- 
tuations which can trigger or enhance shear in another 



location by exceeding a local Coulomb yield criterion. 
A simple probabilistic model of this process describes 
the experimental mean flow across the pipe cross section 
fairly well, but the microscopic picture based on single- 
particle hopping into "holes" seems dubious, in light of 
the present work. ("A particle... will jump to the next 
hole with a probability equal to the probability that the 
stress fluctuation is higher than the threshold" 56].) As 
in the Void Model for drainage, this picture should not 
be taken literally since cooperative motion is likely to 
occur. Perhaps, instead, a spot is triggered to move in 
the opposite direction by stress fluctuations in response 
to other spot displacements. Such interactions are not 
considered in the present model and might be important 
in more general situations. 



D. Forced Shear Flows 

Shear flows can also be driven by moving rough sur- 
faces in granular Coucttc cells 65] . Such force shear flows 
seem quite different from the gravity-driven flows dis- 
cussed above, but some similarities arise, which could be 
signs of the sort of cooperative diffusion described here. 
It is usually found that the shear usually localizes in a 
narrow band at the inner cylinder, as in some pipe flows, 
but recently broad shear ban ds ha ve also been observed 
by Fenistein and van Hecke |l25j ]. If outer part of the 
lower wall rotates with the outer cylinder while the in- 
ner part remains fixed with the inner cylinder, then a 
region of shear smoothly connects stagnant regions on 
either side of the interface. Rising through the Couctte 
cell, the shear region broadens and then drifts toward the 
inner cylinder where it eventually localizes. 

The situation is reminiscent of the broad flow region in 
silo drainage, aside from the direction of the velocity (in 
the horizontal plane of rotation) . In both cases, the veloc- 
ity proflle has a self-similar form reminiscent of diffusion 
problems. In the Couette flow, it is an error function, 
as if connecting two uniform concentrations, rather than 
a Gaussian (for a point source). Unlike gravity-driven 
drainage, the scaling of the shear zone is not that of sim- 
ple diffusion, since its width grows somewhat faster than 
the square root of height, prior to localization near the 
inner cylinder. Nevertheless, it may be that some kind of 
spot-like cooperative diffusion originating with dilation 
at the shearing interface and propagating upward may 
explain the velocity profile of the shear band. 



E. Statistical Thermodynamics of Dense Flows? 

As discussed in the introduction, classical statistical 
thermodynamics does not apply to slow granular flows, 
because the usual definition of temperature, in terms of 
undissipated internal energy, makes no sense. In the op- 
posite limit of static granular systems, some radically 
different definitions have been proposed. For example, 
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Ngan's "mechanical temperature" attempts to describe 
frozen configurational entropy in terms of contact-force 
distributions |25j |. 

Edwards' thermodynamic theory of granular com- 
paction 1^ 1^ F^. I23I focuses instead on excess volume 
and is thus closer in spirit to the Spot Model. His basic 
assumption is that all jammed configurations at a given 
volume are equally probable. This form of ergodicity 
is analogous to Boltzmann's famous postulate that all 
states at a given energy are equally probable. 

Following this analogy, Edwards defines the "com- 
pactivity", analogous to temperature, as the derivative 
of volume with respect to entropy (number of jammed 
configurations). The compactivity can also be related to 
volume fluctuations, e.g. as measured in vibrational ex- 
periments [t^, just as the usual temperature is related to 
velocity fluctuations. Recent granular simulations have 
demonstrated an effective Einstein relation in bulk shear 
flows, where the "temperature", defined as the ratio of 
the diffusivity to the mobility (in response to a force 
pulling on on e g rain), is consistent with the Edwards 
compactivity [63 ■ This suggests that the Edwards en- 
semble may have relevance for dense flows, although it 
is not clear how to predict the diffusivity or derive equa- 
tions for hydrodynamics. 

It is tempting to speculate that the spot density in silo 
drainage is somehow related to the Edwards temperature 
(or compactivity). Indeed, the diffusivity is particles is 
proportional to the spot density, so the Spot Model dis- 
tinguishes between a "warm" region of fast dynamics and 
"cold" region of slow dynamics by the amount of free vol- 
ume. This also seems consistent with viewing the com- 
pactivity as an intensive variable controlling diffusion. 

The drift of particles, however, is also proportional to 
the spot density, and thus, if the analogy were correct, to 
temperature. This troubling implication, which violates 
the notion of "temperature" as a measure of fluctuations, 
is related to the experimental fact that fluctuations de- 
pend on the distance dropped, not time .^j. The reader 
may conclude that any attempt at a thermodynamic de- 
scription of dense granular flows is unlikely to succeed. 

On the other hand, the Spot Model may provide a basis 
for some kind of thermodynamics. Rather than describ- 
ing the system of strongly interacting particles, which 
is wrought with difficulties related to geometrical pack- 
ing constraints, perhaps we should consider a system of 
weakly interacting spots. We have already seen that this 
gives a simple explanation of density waves, but perhaps 
it could also be used to recover thermodynamics. The 
spot diffusivity would be determined by a "spot tempera- 
ture" (perhaps related to the Edwards definition) and the 
spot drift by independent external forcing due to gravity. 
Even if spots were to obey classical statistical mechanics 
(starting from a "spot Hamiltonian" ) , the particles would 
behave rather differently, with drift strongly coupled to 
diffusion. 



F. Compaction Dynamics 

Even in the absence of flow, granular materials exhibit 
nontrivial relaxation phenomena, such as compaction (or 
densification) under vibration. Experiments on vertically 
vibrated granular materials have revealed a very slow 
(logi)~^ decay of free volume with time 75, '126]. To 
explain this result, Boutreux and de Gennes .127.] bor- 
row classical equations from glass theory, due to Cohen 
and TurnbuU l35] (see below): (i) They assume a Poisson 
distribution of "holes" of volume V, 

P{V) = T^e-^/^/ (60) 

where Vf is the local "free volume" , defined (albeit am- 
biguously mill 113) so that Vf = in the hypotheti- 
cal "random close packing" . (ii) The empirical relation, 
Vf oc (logi)~^, then follows from the mean- field assump- 
tion that the rate of compaction, dVf/dt, is proportional 
to the probability 

P{Vn) - e-^°'^' (61) 

of finding a hole larger than a particle volume, Vo. Such 
a hole, capable of being filled by a single particle, corre- 
sponds to a "void" [TSf . 

This thinking reflects the entrenched notion that struc- 
tural rearrangement in an amorphous material occurs 
by the independent displacement of a particle into a 
randomly appearing void. Unlike hypothetical voids in 
granular drainage, which propagate by exchanging places 
with particles, those in granular compaction are annihi- 
lated by hopping particles. This picture of compaction 
has been compared to the "parking-lot model" of cars 
trying to flU a dwindling set of available spaces [t^ ■ 

We have shown that voids do not play a signiflcant role 
in dense granular drainage, so it seems unlikely that they 
could also occur in compaction dynamics, which occurs 
at ever higher volume fractions, closer to the jammed 
state. We would also expect nearest-neighbor cages to 
remain largely preserved during compaction, especially 
in the later stages, when the density is much closer to 
the jammed state than in the drainage flows considered 
here (which exhibit very slow cage breaking). It seems 
more plausible that particles rearrange in a cooperative 
fashion during compaction, slowly eliminating intersti- 
tial free volume without ever opening any particle-sized 
voids. 

The Spot Model provides a general theoretical frame- 
work to describe such collective mechanisms for trans- 
porting free volume from the bulk to the upper surface. 
Note that a new microscopic mechanism need not con- 
tradict the prior successes of mean-field free- volume mod- 
els. As with the mean velocity of granular drainage, the 
mean compaction rate could be very similar with spots 
and voids, even though diffusion and cage-breaking are 
completely different in the two cases. The characteristic 
free volume, wq, in Eq. (|51|) should be related to the typ- 
ical volume carried by a spot, which is roughly a tenth 
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of a particle volume in granular drainage. This is con- 
sistent with the more general view of Vq as the "typical 
free- volume fluctuation required for a local collapse" (see 
below) [13. 



X. OUTLOOK FOR GLASSES 
A. Glassy Relaxation 

A major reason for the current excitement about gran- 
ular materials in physics is the analogy with glasses, 
which exhibit similar features such as "structural tem- 
peratures" (TJj and a "jamming transition" The 
hope is that such common features could s hed l ight on 
the nature of glass and the glass transition Il28| , which 
has been called "the deepest and mos t int eresting un- 
solved problem in solid state theory" |129| |. As in the 
case of granular flow, many toy models and continuum 
approaches have been proposed, but glassy relaxation re- 
mains poorly understood at the microscopic level. 

It has long been recognized that cooperative motion 
plays a crucial role in glassy relaxation, although pre- 
cise mechanisms remain elusive. Early theories of the 
diverg ing r e laxa tion time at the glass-transition temper- 
ature |l30l Il3l| were bas ed on the hypothesis of free- 
volume diffusion m S [m [l3|. Adam and Gibbs 
flrst suggested instead that free volume is associated 
with temperature-dependent regions of cooperati ve re - 
laxation, whose size diverges at the glass transition |l34j . 
Cohen and Grest later introduced more explicit correla- 
tions into free-volume theory (see below) by consider- 
ing clusters of "liquid-like cells", so the melting of glass 
could be u nderstood as a dynamical percolation tran- 
sition |122| . The mode-coupling theory of Gotze and 
Sjogren has invoked the "cage effect" of local blocking 
at the scale of sin gle p articles to explain the emergence 
of slow dynamics 'l35l|, albeit without describing how it 
occurs in real space. In contrast, Berthier has provided 
new simulation evidence for the Adam-Gibbs hypothesis 
that the glass transition is associated with dynamical cor- 
relations at a len gth scale which grows with decreasing 
temperature |l36l |. 

Experiments and simulations have revealed ample 
signs of "dynami cal h eterogeneity" in supercooled liquids 
and spin glasses |l37j , but the direct observation of coop- 
erative motion has been achieved only recently. Rather 
than compact regions of local arrangement, Donati et al. 
have observed "string-like" relaxation in molecula r dy - 
namics simulations of a Lennard- Jones model glass |l38l | . 
The strength and length scal e of correlations increases 
with decreasing temperature |l39j| . consistent with the 
Adam-Gibbs hypothesis. These observations have moti- 
vated Garrahan and Chandler to propose a phenomeno- 
logical lattice model where each site possesses an arrow 
pointing in the "mean direction of facilitation" (opposite 
to the preferred direction of flow) and chain-like relax- 
ation follows from dynamical interactions between the 



arrows |l4Clj |. (The idea of a local vector state variable 
also arises in theories of plasticity, discussed below.) 

Cooperative motion would be difficult to observe ex- 
perimentally in a molecular glass, but Weeks et al. have 
used confocal microscopy to reveal three-dimen siona l 
clusters of fast- moving particles in dense colloids |l4lj| . 
In the supercooled liquid, clusters of cooperative relax- 
ation have widely varying sizes, which grow as the glass 
transition is approached. In the glass phase, the clusters 
are much smaller, on the order of ten particles, and do 
not produce significant rearrangements on experimental 
time scales. 

These observations suggest that spots may have rele- 
vance for cooperative diffusion in simple glasses. String- 
like relaxation is reminiscent of the trail of a spot (see 
Fig. ■ Atomic chains would correspond to the motion 
of spots roughly one particle in size and carrying less 
than one particle of free volume. Larger regions of cor- 
related motion might involve larger spots and/or collec- 
tions of i ntera cting spots. Key features of the experimen- 
tal data |l4lj seem to support this idea: (i) Correlations 
take the form of "neighboring particles moving in paral- 
lel directions" (as in Fig. [7|); and (ii) the large clusters 
of correlated motion tend to be fractals with dimension 
two, as would be expected for the random- walk trail of a 
spot (as in Fig. without the bias of gravity). It seems 
more experiments and simulations should be done to test 
the spot hypothesis for various glass-forming systems. 



B. Free- Volume Theory 

For metallic glasses |l42l Il43j , free- volume theory pro- 
vides the basis for modern continuum models of plastic 
flow and diffusion. The concept of hole diffusion was 
apparently introduced in 1936 by Eyring |9^, but Co- 
hen and Turnbull in 1 959 initia ted the standard model 
of glassy relaxation [95lll44 fl45j . based on Eqs. ll^ and 
(|61|1 . In their theory of "molecular transport in liquids 
and glasses" (on the same footing), randomly distributed 
"holes" are filled by single-particle displacements, when- 
ever a hole is large enough to accommodate a neighbor- 
ing particle, with probability P{vq) — e^°^^f . These as- 
sumptions yield the classical formula, D (x e~^°/'"f ^ for 
the diffusivity, where a linear dependence of free volume 
on temperature, Vf = A + B{T — Tg), above the glass 
transition, Tg, is commonly assum ed, as suggested by 
Williams, Landel, and Ferry |l3lj |. By assuming that 
hole diffusion is biased by the local shear stress, Spaepen 
derived Doolittle's formula for the viscosity, rj oc e'"°^^f , 
and proposed mechanisms for the creation and annihila- 
tion of free volume based on particles jumping into holes 
of a smaller or larger volume, respectively 146]. 

Of course, this microscopic picture is very similar to 
that of the Void Model for granular drainage, where grav- 
ity provides a bias for void diffusion (so voids are created 
at the orifice and de stroy ed at the free surface). Indeed, 
Spaepen's Figure 3 |l46l | is nearly identical to Figure [21 
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above. The key point is that, in both cases, flow and 
diffusion occur when a single particle moves into a void 
independently from its neighbors. 

Since cooperative diffusion is widely believed to oc- 
cur in glasses, the Spot Model seems to provide a more 
realistic mechanism for free volume diffusion. As we 
have demonstrated for granular drainage, this could be 
checked in glasses by comparing the diffusion lengths for 
free volume and particles and the cage-breaking length 
(which must all be roughly the same for hole diffusion) or 
by directly observing spatial velocity correlations. If in- 
teracting spots exist in metallic glasses, perhaps they are 
responsible for the formation of shear bands of reduced 
density under large tensile loading |147| . Shear bands in 
metallic glasses are difficult to understand, because, un- 
like shear bands in crystals, which are known to result 
from dislocation interactions, it is not clear what are the 
microscopic carriers of plastic deformation in disordered 
materials. 



C. Plastic Deformation 

The idea of localized cooperative motion has a long 
history i n th e plastic deformation of amorphous solids. 
Orowan |l48j | was perhaps the first to propose, contrary 
to Erying's hole-diffusion hypothesis, that plastic fiow 
occurs via localized shear transformations in regio ns of 
enhanced atomic disorder (and free volume) |148| |. In 
the context of metallic glasses. Argon Il49| later distin- 
guished between two types of local rearrangements: (i) 
a "diffuse shear transformation" at high temperature, in 
which a spherical region 4-6 atoms in diameter deforms 
smoothly into an ellipsoidal shape, and (ii) an "intense 
shear transformation" at low temperature, in which an 
equally small disk-shaped region, similar to a disloca- 
tion loop, slides suddenly along an atomic slip plane. 
A generalization of this concept, the "localized inelastic 
transformation" (LIT), forms the basis of a stochastic 
model of ela sto-plastic deformation due to Bulatov and 
Argon |l5Clj . Each LIT acts irreversibly on a material 
element as a pure shear along some discrete axis of sym- 
metry plus an additional dilation. 

The basic ide a of the LIT has recently been re- 
vived i n ph ysics |l5lj as a "shear transformation zone" 
(STZ) '152], connecting "-h" and "-" states relative to 
the direction of an applied shear. Falk and Langer first 
identified STZ events in molecular dynamics simulations 
of a Lennard-Jones glass and proposed a new continuum 
theory of viscoplastic deformation, which involves the 
densities, n± , of the two hypothetical STZ types as aux- 
iliary continuum variables (or order parameters) |l52j| . 
They noted that STZ transition rates should depend on 
the local free volume, taken to be constant as a first ap- 
proximation. STZ theory compar es fav orably with exper- 
imental data on metallic glasses |l53j . although it con- 
tains a number of adjustable parameters and ad hoc func- 
tional forms. 



Lemaitre has recently extended the two-state STZ the- 
ory to account for the creation and annihilation (but not 
diffusion) of free volume 60] , borrowing some equations 
from the literature on metallic glasses and granular com- 
paction discussed above. The resulting empirical consti- 
tutive laws are meant to apply to a broad class of sheared 
dense materials, including glasses [154| and granular ma- 
terials _61j. The same phenomenology provides a fairly 
simple way of understanding universal features of amor- 
phous systems, such as stick-slip, dilatency, and nonlin- 
ear rheology, although again the microscopic mechanisms 
remain vague. Moreover, if free volume is created and 
destroyed by particle displacements, it seems it should 
also di ffuse by very similar mechanisms (as in Spaepen's 
theory [Ti^ '). 

Although neglecting diffusion, Lemaitre offers a new, 
qualitative interpretation of free- volume kinetics, of par- 
ticular interest here: "Borrowing on the argument of 
STZ theory, compaction results from elementary rear- 
rangements, involving several molecules at the meso- 
scopic s cale" , rather than "the motion of a single grain in 
a hole" |l54j . Of course, we have reached the same con- 
clusion in the completely different context of granular 
flow. Lemaitre gives no further description of the "el- 
ementary rearrangements" associated with free-volume 
dynamics, but the Spot Model provides a possible math- 
ematical framework, which also naturally accounts for 
diffusion. This seems like an interesting direction for fu- 
ture research. 



XI. CONCLUSION 

Returning to our original question about random 
walks, we arrive at a simple dynamical distinction be- 
tween different states of matter: 

Gas. The particles in a gas undergo independent ran- 
dom walks, as ballistic trajectories are randomly 
redirected by collisions. 

CrystaL The particles in a (single) crystal diffuse by 
hopping to either interstitial or vacant lattice sites. 
Both interstitials and vacancies undergo thermally 
activated, independent random walks on the lattice 
(at sufficiently low defect density). 

Liquid. The particles in a liquid move cooperatively at 
the scale of several particles, but, due to the high 
internal kinetic energy, they undergo independent 
random walks over much longer distances at exper- 
imental time scales. 

Amorphous. The particles in an amorphous state move 
cooperatively at experimental time scales, due to 
the low internal energy. Proximity to jamming pre- 
cludes the vacancy /interstitial mechanism (insert- 
ing particles into holes). Instead, particles undergo 
locally correlated random walks in response to dif- 
fusing spots of free volume. 
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This is surely an oversimplificatioii. For example, direc- 
tional covalent bonding can significantly alter th e pict ure: 
In amorphous silicon, vacancy (dangling bond |l5Rl |) or 
interstitial (floating bond 15fr|) mechanisms may domi- 
nate relaxation, while in open lattices like diamond sili- 
con, coop erati ve self-difFusion can also occur (concerted 
exchange |l57j ). Nevertheless, for materials with dense 
disordered phases, the state classification above by diffu- 
sion mechanism seems reasonable. 

In summary, the Spot Model provides a simple math- 
ematical framework for cooperative diffusion, when in- 
dependent random walks are inhibited by packing con- 
straints. Like the original random- walk concept, which 
is now used in many fields other than physics, the ba- 
sic idea of a moving "spot of influence" may also have 



broader applications. 
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